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An efficient continuous-time path-integral Quantum Monte Carlo algorithm for the lattice polaron 
is presented. It is based on Feynman's integration of phonons and subsequent simulation of the 
resulting single-particle self-interacting system. The method is free from the finite-size and finite- 
time-step errors and works in any dimensionality and for any range of electron-phonon interaction. 
The ground-state energy and effective mass of the polaron are calculated for several models. The 
polaron spectrum can be measured directly by Monte Carlo, which is of general interest. 
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The past few years have witnessed a rapid develop- 
ment of continuous- time Quantum Monte Carlo (QMC) 
algorithms for quantum- mechanical lattice models. The 
driving force behind this is the desire to eliminate the sys- 
tematic errors introduced by the Trotter decomposition 
in the standard discrete-time QMC methods |Q. The 
main idea then is to regard the imaginary-time evolu- 
tion of a particle or spin configuration as a continuous- 
time Poisson process with "events" being either a particle 
jump or a spin flip. In this way continuous-time QMC 
algorithms were developed for a particle in an external 
potential Heisenberg model HQ, t ~ J model M, 
bosonic Hubbard model B, and Frohlich polaron [Q. 

In this Letter I present a continuous-time path-integral 
QMC algorithm for the lattice polaron, i.e., an elec- 
tron strongly interacting with phonons on a lattice. The 
method combines analytical integration of the phonon de- 
grees of freedom with the advantages of the continuous- 
time formulation of the Monte Carlo process. The 
method is universal. It works for infinite lattices in any 
dimensionality and for any radius of electron-phonon in- 
teraction. It is also free from the systematic finite-time- 
step errors which were an undesirable feature of the orig- 
inal (discrete-time) QMC algorithm based on the inte- 
gration of phonons In fact, the new method al- 
lows for exact (in the QMC sense) calculation of the 
ground-state energy, effective mass, and even spectrum 
of the polaron. The numerical accuracy of the method is 
0.1 — 0.3%, which is not as good as that of exact diagonal- 
ization PJlfj| and density-matrix renormalization group 
p| schemes, but is good enough for practical purposes. 

I begin with the important question of which quanti- 
ties can be calculated with a path-integral QMC method. 
The traditional scheme samples paths which are peri- 
odic in imaginary time, see, e.g., Refs. jjlp|. This allows 
for the calculation of thermodynamic properties such as 
internal energy or specific heat, as well as some static 
correlation functions, but not dynamic properties of the 



system. In Ref. |12[ it was shown that the polaron effec- 
tive mass, an important dynamic characteristic, can be 
measured on the ensemble of paths with open boundary 
conditions in imaginary time (BCIT). Here I extend the 
argument to the whole polaron spectrum. Consider two 
quantities. The first one, Zp = J2 i {i\e~ l3H \i)5p^p, is 
by definition the partition function of many-body states 
with fixed total (electron plus phonon) quasimomentum 
P [(3 = (/csT)" 1 and H is the Hamiltonian] . The sec- 
ond quantity is the partition function with twisted BCIT, 
Zat = TrA r e _/3ff , in which the many-body configura- 
tions (electron position r and ionic displacements £) at 
t — (3 are obtained from the configurations at t = by 
shifting along the lattice by the vector Ar. There exists a 
fundamental Fourier-type relation between the two fl^| : 



Ar 



<z 



Ar 



VvVie^^wW)i{r)\. (1) 



Here w[r(r)^(r)] is the (positive-definite) weight of 
a many-body path satisfying the twisted BCIT, and 
J VrV^ means the integration over such paths with 
alf possible shifts Ar. In the low-temperature limit, 
f3 — > oo, Zp is dominated by the lowest eigenstate Ep: 
Z P =exp(-/3£p) and Ep=-^?§£-. Equation (0) then 
gives 
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where (A) = {J^VrV^w)- 1 J op VrV^Aw stands for 
the average in the case P = 0. Equation (||) shows that 
the ground-state dispersion of the polaron can be inferred 
directly from imaginary-time simulations. However, the 
formalism is not free from the sign-problem, as apparent 
from Eq. (Q). Energies for non-zero values of P can be 
obtained only at intermediate and large electron-phonon 
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couplings when the average (e lPAr )o = (cosPAr)o is not 
small. The ground state corresponds to P = and its 
energy can be calculated straightforwardly using Eq. (|^). 

Analogous formulas can be obtained for derivatives of 
Ep with respect to momentum, e.g., for the inverse ef- 
fective mass. In the low-temperature limit one can write 
Ep = — 4 In Zp since there is no difference between 1/(3 
and d/d/3. From this and Eq. (Q), it follows that 
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One can see that the inverse effective mass is the "diffu- 
sion coefficient" of the imaginary-time propagation. (For 
similar treatment of the continuous polaron and of an 
atom of 3 He in liquid 4 He, see |l3| and H], respectively.) 

Thus, to compute dynamical properties, the QMC pro- 
cess should sample paths with open BCIT. Consider now 
a particular polaron system on a hypercubic lattice with 
nearest-neighbor hopping. The model is defined by the 
Hamiltonian 
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The phonon subsystem (operators b^, b m ) is a set of un- 
coupled harmonic oscillators, one per site, with inter- 
nal coordinates £ m , frequency to, and reduced mass M. 
The lattice is assumed to be infinite in all dimensions. 
The electron-phonon interaction is taken to be of the 
"density-displacement" form. No restriction is imposed 
on the form of force / m (n) with which the electron at 
site n acts on m'th oscillator. 

Since Feynman's classic work on polarons [ jDs] , it has 
been known that phonon degrees of freedom (variables 
£(t)) can be integrated out in the path-integral repre- 
sentation of the partition function, Eq. (|l]), leading to 
a retarded self-interaction of the electron. For periodic 
BCIT, £ m (/3) = £m(0), the phonon-induced part of the 
polaron action is |15| 
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where r(r) is the electron path. However, in the case of 
open BCIT the integration must be performed under the 
constraint £ m +Ar(/3) = £m(0), where Ar = r(/3) — r(0) is 
the shift of the electron path. For nonzero Ar this yields 
an extra term in the action 



AA[r(r)] = — £s m (C m+Ar - C n 
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where 
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This formula is valid only in the limit e 13 3> 1, which 
is easy to realize in practice. Thus, the integration over 
phonons reduces the problem to a single-particle system 
with an extra factor exp(A) — exp(A per + AA) in the 
weight of each path r(r). 

I now describe how this factor should be incorporated 
into the general scheme of continuous-time QMC U& . In 
d dimensions the hopping term in the Hamiltonian (H) in- 
troduces 2d independent Poisson processes with events, or 
"kinks" (after Ref. M), being jumps of an electron path 
to one of 2d nearest neighbors. The probability to have 
N k kinks of a given sort on the time interval [0, (3\ is given 
by the Poisson distribution P Nk = {N k \y 1 {tf3) Nk e~ . 
The Monte Carlo process consists of (i) proposing a 
change of the path by either removing existing kinks or 
adding new kinks and (ii) either accepting or rejecting 
the proposal. (Another possible subprocess is the shift of 
a kink in time but this can always be achieved by adding 
and removing.) Because of the open BCIT, it is sufficient 
to change the number of kinks just by one. When a kink 
is added (removed) at time To the whole path at t > To 
is shifted in the corresponding direction (antidirection) 
by one lattice site (see Fig. |j]) . The balance equation for 
the adding-removing process is 



q a W(N k )P a (N k -> N k + 1) = 



q r W(N k + l)P r (N k + 1 -» N k ), 
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where W(N k ) is the probability to have a given path with 
N k kinks of a given sort, and q a and q r are the probabil- 
ities of attempting to add or remove a kink, respectively. 
In this paper, q a =q r — 1/2 is used for N k >l, and q a =l, 
q r =0 for N k =0 (if there are no kinks of a given sort, one 
can only add one). In the absence of electron-phonon 
interaction the ratio W(N k + l)/W(N k ) = t/3/{N k + 1) 
follows from the Poisson distribution. In the general case, 
each W is multiplied by its phonon-induced weight e A . 
The acceptance rules now follow from Eq. (@): 
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to add a new kink to the existing N k of a given sort, and 

A fe + 1 



Pr{N k 



1 



N k ) = min 



tj3 



„A Nk -A Nk + 1 



(10) 



to remove one of the existing + 1 kinks of a given sort. 
For the special case N k — 0, the preexponential factor in 
Eq. (g) must be t(3/2 instead of tfi, and in Eq. |l(]) it 
must be 2/(t/3) instead of l/{t(3). 
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FIG. 1. A typical polaron path for /3 = 10t _1 (solid line). 
When the kink at to = 7.0 1" 1 (bold line) is removed, the 
whole path at to < r < (5 is shifted to the left by one lattice 
site (dashed line). Ar changes from -1 to -2. 



The explicit form of the energy estimator that enters 
Eq. (||) follows from the Ar — > limit of the correspond- 
ing finite-time expression S 
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where N^ ot is the total number of kinks (of all sorts) on 
a path. Note that due to the open BCIT the thermody- 
namic estimator ( pd| ) measures the ground-state energy 
rather than full internal energy. 

Summarising the procedure, the QMC process sam- 
ples single-particle paths with open BCIT by inserting 
and deleting kinks of 2d types. The acceptance rules for 
these two fundamental processes are given by Eqs. (||) 
and (|lo|). The action A is a functional of the electron 
path r(r) and is given by the sum of Eqs. ([|) and 
Measured quantities include the energy estimator ft\A\ ) 
and the inverse mass estimator 4(r a (/3)— r Q (0)) 2 . The 
polaron spectrum and effective mass are calculated with 
Eqs. (||) and (||), respectively. Note that statistics for all 
momenta (i.e., for as many P-points as one wants) are 
collected during a single QMC run. In practical simu- 
lations successive measurements were taken every tenth 
single-kink step to reduce statistical correlations. For 
each set of model parameters several series of 500 000 or 
1 000 000 measurements were conducted for different val- 
ues of j3 to detect possible finite-temperature systematic 
errors, typically for (3Huj — 10, 15, 20, and 25. No such 
errors were detected. 

The method was tested on the Holstein model, for 
which independent numerical results are available. The 
Holstein model is a particular case of Eq. (|J) with 
/m(n)=ft<5 mn . The model is parametrized by the 
dimensionless frequency uj=hu)/t and the dimension- 
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FIG. 2. Ground-state polaron energy for several models. 
In all cases, Q — 1.0. Statistical errors are smaller than the 
symbols. 



less coupling constant A = (J2 m f^G)) / (AdtMuj 2 ) = 
re 2 / (AdtMuj 2 ). Excellent agreement with previously pub- 
lished or communicated results was found. For example, 
in d = 1 [Figs. || and || (circles)] the present method's 
estimate of the ground state energy for Q = 1.0, A = 0.5 
is E = -2.471 ± 0.001 (in units of t), with the true 
value being in the interval (-2.46968, -2.471) @; for 
Q = 1.0, A = 1.0 it yields E = -2.999 ± 0.001 with 
the true value in (-2.99883, -3.000) @. For Q = 1.0 and 
A = 1.25, 1.5, and 2.0, the estimated energies are -3.298, 
-3.623, and -4.388, respectively, which are, to within sta- 
tistical error (±0.002), precisely the values obtained by 
the exact diagonalization (ED) method Jig ]. For Cj = 2.0 
and A = 1.5625 and 2.25, QMC yields E = -4.013 
and -5.070, respectively, in agreement with the strong- 
coupling perturbation theory p9| . Polaron masses ob- 
tained by QMC arc in excellent agreement with varia- 
tional calculations |l^1 and with density-matrix renor- 
malization group (DMRG) results jn|. In d = 2, ui = 1.0 
[Figsj2 and || (filled squares)], QMC's energies agree with 
ED |g and QMC's masses agree with DMRG @. All 
of these checks confirm that the present algorithm is free 
from systematic errors of any kind. Statistical errors are 
small, 0.1% — 0.3% in most cases. 

As a new application of the method the ground-state 
energy and effective polaron mass were calculated for the 
three-dimensional Holstein model [see Figs. | and | (t ri- 
angles)]. At present, three-dimensional lattices are out of 
reach of ED and DMRG methods due to the enormous 
phonon Hilbert space. The present Monte Carlo algo- 
rithm treats all the phonons in an infinite lattice exactly 
through the analytical integration so the lattice size is 
irrelevant. Another advantage of the method is its abil- 
ity to study long-range electron-phonon interactions. In- 
deed, the interaction enters the formalism only via lattice 
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FIG. 3. Inverse effective polaron mass in units of 
1/mn = 2ta 2 /Tl 2 for several models. In all cases, Q = 1.0. 
Statistical errors are smaller than the symbols. 
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FIG. 4. The polaron spectrum Ep (in units of t) in the 
one-dimensional Holstein model for uj = 1.0 and A = 1.75. 
Open circles are exact diagonalization results [20]. 



sums ^ m /m(i"i)/m(r2) which can be computed in ad- 
vance for all ri — T2 and stored for later calculation of the 
action A. Figures || and || (open squares) show results for 
im(n) = «(|m— n| 2 + l)~ 3 / 2 in d = 2 (for this interaction 
A = 1.742 K 2 /(8tMu> 2 )). This form of / m (n) has recently 
been proposed to model the interaction between in-plane 
holes and apical oxygens in high-T c superconductors pl| . 
One can see that at large A the polaron is much lighter 
than in the Holstein case. Indeed, a long-range interac- 
tion results in "wider" polaron paths and, consequently, 
in a smaller effective mass by virtue of Eq. (|J). Finally, I 
present the polaron spectrum in the d = 1 Holstein model 
for Q = 1.0 and A = 1.75, calculated with Eq. (||) (see 
Fig. |^) . The spectrum flattens at large momenta as was 
revealed in previous numerical studies [^2|,^3| . The actual 
energies Ep are in excellent agreement with ED results 
pof . This demonstrates that a real-time spectrum can 
be measured directly by an imaginary-time path-integral 
QMC. 

In conclusion, an efficient continuous-time path-integ- 
ral Quantum Monte Carlo algorithm for the lattice po- 
laron has been presented. It is free from the finite- 
size and finite-time-step systematic errors. The method 
works for infinite lattices and any range of electron- 
phonon interaction. It allows for exact calculation of the 
ground-state energy, effective mass and spectrum of the 
polaron. More technical details will be presented else- 
where. 
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